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Abstract 

Background: In this paper we propose a model reduction method for biochemical reaction networks governed by a 
variety of reversible and irreversible enzyme kinetic rate laws, including reversible Michaelis-Menten and Hill kinetics. 
The method proceeds by a stepwise reduction in the number of complexes, defined as the left and right-hand sides 
of the reactions in the network. It is based on the Kron reduction of the weighted Laplacian matrix, which describes 
the graph structure of the complexes and reactions in the network. It does not rely on prior knowledge of the 
dynamic behaviour of the network and hence can be automated, as we demonstrate. The reduced network has fewer 
complexes, reactions, variables and parameters as compared to the original network, and yet the behaviour of a 
preselected set of significant metabolites in the reduced network resembles that of the original network. Moreover 
the reduced network largely retains the structure and kinetics of the original model. 

Results: We apply our method to a yeast glycolysis model and a rat liver fatty acid beta-oxidation model. When the 
number of state variables in the yeast model is reduced from 1 2 to 7, the difference between metabolite 
concentrations in the reduced and the full model, averaged over time and species, is only 8%. Likewise, when the 
number of state variables in the rat-liver beta-oxidation model is reduced from 42 to 29, the difference between the 
reduced model and the full model is 7.5%. 

Conclusions: The method has improved our understanding of the dynamics of the two networks. We found that, 
contrary to the general disposition, the first few metabolites which were deleted from the network during our 
stepwise reduction approach, are not those with the shortest convergence times. It shows that our reduction 
approach performs differently from other approaches that are based on time-scale separation. The method can be 
used to facilitate fitting of the parameters or to embed a detailed model of interest in a more coarse-grained yet 
realistic environment. 

Keywords: Kinetic models, Enzyme kinetics, Complex graph, Weighted Laplacian, Yeast glycolysis, Rat liver beta 
oxidation 



Background 

A kinetic model of a biochemical reaction network con- 
sists of a set of ordinary differential equations describing 
the dynamics of the concentrations of all metabolites in 
the reaction network. Most biochemical reaction net- 
works are complex and involve many enzyme-catalyzed 
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processes with non-linear kinetics and intricate stoichio- 
metric and regulatory interactions between the enzymes. 
Consequently, the mathematical models of such networks 
contain high-dimensional sets of coupled rational dif- 
ferential equations, which sometimes require huge com- 
putational effort to analyze. The current state-of-the-art 
numerical tools for stability analysis, for bifurcation study 
and for other types of dynamical analysis are known to 
suffer from a so-called curse- of- dimensionality. For exam- 
ple, the largest biological model that has numerically been 
analyzed for bifurcation in [1] consists of 25 metabolites 
and 37 parameters and the one in [2] has 22 metabolites. 
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Moreover since complex models of biochemical reaction 
networks involve a huge number of parameters, the task of 
identifying these parameters (in addition to those parame- 
ters that have been identified in the literature) is enormous 
and requires large datasets. The complexity of this task 
is further compounded by the fact that often not all the 
metabolite concentrations can be measured. Thus, there 
is a need for techniques that can reduce a given kinetic 
model of a biochemical reaction network to a simplified 
version that mimics the behaviour of the original model 
satisfactorily, but contains less differential equations and 
parameters. 

For biochemical reaction networks, a number of model 
reduction techniques are known. See [3] for a detailed 
review of some of the well known methods of model 
reduction. Here we list only a few of the known 
methods in the literature. The singular perturbation 
method, the time-scale separation technique [4-8], the 
rapid-equilibrium approximation, also known as the 
quasi-equilibrium approximation (see [9]) and the quasi 
steady- state approximation (QSSA) (see for e.g., [10]) are 
the most commonly used techniques. The reduced state 
vector obtained by any of these techniques contains a 
subset of the metabolite concentrations (state vector com- 
ponents) of the full model. Hardin [11] extends the QSSA 
approach by considering the higher-order approximation 
in the computation of the quasi steady-state. In [12], some 
approaches for reduction of complexity in biochemical 
reaction networks are considered for use in SYCAMORE, 
which is a computational research environment in sys- 
tems biology. One of the approaches considered in [12] is 
the intrinsic low-dimensional manifold (ILDM) approach, 
which is an improved time-scale separation technique. 
More recently [13], a computational singular perturbation 
(CSP) algorithm was developed to analyze the time-scales 
of the NF-/cB signaling network which could be used in 
order to reduce its model. The application of singular per- 
turbation, time-scale separation, quasi equilibrium and 
quasi steady state approaches to general enzyme-kinetic 
rate laws, such as Michaelis-Menten and ping-pong bi-bi 
is difficult and leads to complicated rate law expressions 
in the reduced models. Some of the time-scale separation 
techniques are either based on a priori experimental infor- 
mation or eigenvalues of the Jacobian corresponding to 
the model and hence are only locally effective. Another 
recent approach for model reduction uses tropical geom- 
etry (see e.g., [14]), wherein the polynomial occuring in 
every rate equation is replaced by a monomial which is 
equal to the largest, in absolute value among the monomi- 
als that constitute the polynomial. 

One of the ways of reducing the complexity of a bio- 
chemical model is to reduce the number of parameters 
in it. This can be done by carrying out a parameter sen- 
sitivity analysis (see for e.g., [15]) and eliminating those 



parameters whose variations have least effect on the 
dynamics of a given network. In [8], a time-scale analysis 
is first done using experimental data to identify the fast 
and the slow metabolites and this information is then used 
to carry out an a priori parameter sensitivity analysis to 
obtain a reduced kinetic model of a biochemical reaction 
network. In [16], an implicit multiparametric variabil- 
ity analysis (MPVA) method is used to search the entire 
parameter space in order to determine the set of param- 
eters that can be eliminated. In [17], the authors go one 
step further by identifying a region in the parameter space 
where some of the parameters are zero-valued, others 
have readjusted values and where nevertheless the outputs 
of the original model match those of the reduced model 
within a certain tolerance. They further show that param- 
eter sensitivity analysis approach for model reduction may 
not be always successful. Another way of reducing the 
complexity of a biochemical reaction network model is to 
reduce the number of reactions. In [18-20], optimization 
techniques are used to determine which reactions can be 
eliminated from the original model without a substantial 
alteration of the model behaviour. 

The method proposed in [21] simplifies a given chem- 
ical reaction network by linearly combining reactions in 
such a way that the resulting network has lesser number 
of species. However, the kinetics for the reduced set of 
reactions are determined by the rate limiting steps in the 
original network and this requires prior biological knowl- 
edge of the network. Dano et al. [22] propose reduction by 
identification and elimination of variables in such a way 
that the basic dynamic properties of the original model 
are preserved. In [23], model reduction is achieved by 
simplifying the rate equations of individual enzymes in 
the network. A limitation of this method is that it yields 
a reduced model that predicts measurement data only 
under specific experimental conditions. 

In this paper, we describe a new model reduction 
method that reduces the number of reactions, metabolites 
and parameters, such that the dynamics of the metabo- 
lite concentrations of the reduced model are close to those 
of the original model. This method proceeds by a simple 
stepwise reduction in the number of complexes', which 
are defined as the left and right-hand sides of the reac- 
tions in the network. The effect of this stepwise reduction 
is monitored by an error integral, which quantifies how 
much the behaviour of the reduced model deviates from 
the original. The method is based on the reduction of 
the underlying weighted Laplacian (see [24] for a defi- 
nition) describing the graph structure of the complexes 
and reactions in the chemical reaction network. A similar 
technique is also employed in the Kron reduction method 
for reduction of resistive electrical network models [25] . 

Our model reduction technique is easy to implement 
and can be used to reduce reaction networks governed by 
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a vast majority of enzyme kinetic rate laws including Hill 
kinetics and reversible Michaelis-Menten kinetics. It does 
not rely on prior knowledge about the dynamic behaviour 
or biological function of the network. Consequently, it can 
be automated. Furthermore, the reduced model largely 
retains the kinetics and structure of the original model. 
This enables a direct biochemical interpretation and yields 
insight into which parts of the network have the highest 
influence on its behaviour. It also accelerates computa- 
tions and facilitates parameter fitting, especially when we 
deal with models of huge biochemical reaction networks. 

We show the application of our model reduction tech- 
nique to a yeast glycolysis model [26] and a model of 
beta-oxidation in rat liver [27]. We have simulated the 
transient behaviour of the metabolites that are not elim- 
inated during the model reduction procedure. In both 
the cases, a 30% reduction of the number of variables 
still retained about 92.5-96.5% agreement between the 
outputs of the full and the reduced networks. 

Methods 

Preliminaries 

Notation: The space of n dimensional real vectors is 
denoted by R n , and the space of m x n real matrices 
by R mxn . The space of n dimensional real vectors con- 
sisting of all strictly positive entries is denoted by R" . 
l n denotes an identity matrix of dimension n. The trans- 
pose of a matrix A is denoted by A T . Define the mapping 
Ln : R™ R m , x i-> Ln(#), as the mapping for which 
the i-th component is given by (Ln(x))i := ln (x£). Simi- 
larly, define the mapping Exp : R m -> R+, x i-> Exp(^), 
as the mapping for which the /-th component is given 
by (Exp := exp (pa). t m denotes a vector of dimen- 
sion m with all entries equal to 1 and dim(l/) denotes the 
dimension of a set V. 

Chemical reaction network structure 

In this section, we introduce the concept of a complex 
graph which was first introduced in the work of Horn 
& Jackson and Feinberg [28-30]. This concept will be 
used first in deriving our general mathematical formula- 
tion of the dynamics of chemical reaction networks, and 
subsequently to explain our model reduction approach. 

Let m, c and r denote the number of species (metabo- 
lites), complexes and reactions respectively of a given 
chemical reaction network. The set of complexes of a 
network is simply defined as the union of all the differ- 
ent left- and righthand sides (substrates and products) of 
the reactions in the network. Thus, the complexes corre- 
sponding to the network in Figure 1 are 2X\ + X2, X3, 
Xi+2X 2 andX 4 . 

The expression of the complexes in terms of the chemi- 
cal species is formalized by an m x c matrix Z, whose a-th 
column captures the expression of the a-th complex in the 



m chemical species. For example, for the network depicted 
in Figure 1, 

"2 0 10" 

10 2 0 

0 10 0 ' 
_0 0 0 1 _ 

The matrix Z is called the complex stoichiometric matrix 
of the network. Note that by definition all elements of the 
matrix Z are non-negative integers. 

Since complexes are the left- and righthand sides of 
reactions in the network, they can be naturally associated 
with the vertices of a directed graph Q with edges corre- 
sponding to the reactions. Formally, the reaction a — > 
between the a-th (reactant) and the /3-th (product) com- 
plexes defines a directed edge with tail vertex being the 
a-th complex and head vertex being the /3-th complex. 
The resulting graph will be called the complex graph. 

Any graph is defined by its incidence matrix B [24] . This 
is a c x r matrix, where c denotes the number of vertices, 
r denotes the number of edges and the (a,;)-th element is 
equal to — 1 if vertex a is the tail vertex of edge j and 1 if 
vertex a is the head vertex of edge while 0 otherwise. 

The basic structure underlying the dynamics of the vec- 
tor x e R+ of concentrations x» i — 1, • • • , m, of the 
species of a chemical reaction network is given by the bal- 
ance laws x = ZBv(x); the elements of the vector v e R r 
are commonly called the (reaction) rates or fluxes. 

In many cases of interest, especially in biochemical reac- 
tion networks, chemical reaction networks are intrinsi- 
cally open, in the sense that there is a continuous exchange 
with the environment (in particular, inflow and outflow 
of chemical species and connection to other reaction net- 
works). In this paper, we are particularly interested in 
inflows to and outflows from the complexes of the net- 
work. This will be modeled by a vector Vb(x) £ R c 
consisting of c boundary (or, exchange) fluxes, leading to 
an extended model 

x = ZBv(x) + Zvb(x) (1) 

A linkage class of a chemical reaction network is a maxi- 
mal set of complexes {Ci, . . . , C/J such that C; is connected 
by reactions to Cj for every i,j e {1, . . . , k}, i 7^ j. 

General kinetics 

For a biochemical reaction network, the relation between 
the reaction rates and species concentrations depends 
on the mechanisms of the reactions involved in the net- 
work. In this section, we derive a framework for describing 
the dynamics of enzymatic reaction networks using the 
aforementioned notion of complex graphs in a reaction 
network. This framework will be useful in describing our 
model reduction method. We describe the forward and 
reverse rates of an enzyme with separate equations. 
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2X 1 + X 2 



Figure 1 Example of a reaction network. 



X 



X, 



X x + 2X 2 



and equation (4) can be written 



Let Zsj denote the column of the complex stoichiomet- observe that d : 
ric matrix Z corresponding to the substrate complex Sj as 
of the y-th reaction of a chemical reaction network. Then 

the unidirectional, forward reaction rate of the ; th reaction v(x) = d(x)kxix 2 = d(x)k exp (zjLn(#)). 

of the chemical reaction network between the substrate 
complex Sj and the product complex Vj is given by 



Vj{x) = dj(x)kjexp (Z^.Ln(x)) t 



(2) 



where for; = l,...,r,dj : —> R + is a rational function 
of its argument and kj denotes a proportionality constant 
of the / h reaction of the network. Note that if the gov- 
erning law of the ; th reaction is mass action kinetics, then 
dj(x) — 1. 
As an example, consider the reaction 



X\ + X2 — > X3 + X4 



(3) 



For i = 1, . . . , 4, let X{ denote the concentration of the 
species X{ and define x :=[x\ X2 X3 X4] 1 '. In this 
example, the substrate complex S is X\ +X2 and the prod- 
uct complex P is X3 + X4. The matrices B and Z for the 
reaction (3) are given by 

"10' 
1 0 
0 1 
0 1 ^ 

Hence Zs is given by Zs = [1 1 0 0] r . Let K\> I<2, 
K3 and 7<4 denote the "Michaelis" constants of the species 
X\, X2, X3 and X4 respectively. Let Vf denote the maxi- 
mum rate of the forward reaction (3). The expression for 
the rate of the reaction (3) depends on the type and mech- 
anism of the reaction, for instance the type of inhibition 
and other regulatory mechanisms involved in the reaction. 
One possibility is 



v(x) = 



(Ah 

f 1 I Xi 1 X3 \ ( 1 - X2 I X4 \ 

\ L Ki K 3 A K 2 ^ K 4 ) 
By defining k := 

d(x) := ? ^- v , 

\ l + + iC 3 A iC 2 + IQ ) 



(4) 



In Additional file 1, we have provided a list of enzyme- 
kinetic rate laws that can be written in the form (2) with d 
satisfying d : R+. Note that the form (2) is retained 

with a d satisfying d : R^ -> R + even if there are compet- 
itive, non-competitive or uncompetitive modifiers in the 
reaction network. This is because the terms containing the 
concentration of the modifiers appear in the denominator 
of the rate expression in such a way that the denominator 
assumes positive values for positive values of concentra- 
tions of the substrates. For example, if we consider the 
reaction 



x 2 



(5) 



governed by Michaelis -Menten kinetics, with Michaelis 
constant of X\ denoted by K\ and the maximum reaction 
rate denoted by V, then 



v(x) 



Vxi 



1 , xi 
1 + Ki 



(6) 



is the expression for the rate of the reaction (5). It is easy 

to see that with k := ^ and d(x) := (l + , (6) has 

the same form as (2) and d : R+ -> R+. Now assume that 
the reaction involves a competitive modifier / whose con- 
centration is denoted by i and whose inhibition coefficient 
is denoted by 7<Q. Then the rate of the reaction is given by 



v c (x) 



Vxi 
Ki 



I 1 xi 
1 ^ Ki 



(7) 



Ki 



With k defined as earlier and d(x) := (l + ^ + ±\ 

observe that (7) has the same form as (2) and d : R+ 
R+. Similarly if / is a non-competitive modifier, then 



Vxi 
Ki 
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denotes the rate of the reaction which can again be writ- 
ten in the form (2) with a d satisfying d : R+. If / 
denotes an uncompetitive modifier, then 



= (*) 



V5ci 
Ki 



i . *L 
1 ^ iCi 



denotes the rate of the reaction which can be written in 
the form (2) with a d satisfying d : R+ R+. 

Based on the formulation in (1) and (2), we can describe 
the complete reaction network dynamics as follows. For 
every a, n G {1, . . . , c}, define 

Cna :={je{l,...,r}\(a,7t) = (Sj,Vj)} 

and ana := ^2j G c na kjdj(x). Thus if there is no reaction 
a — »> 7r, then = 0. Define the weighted adjacency 
matrix A of the complex graph as the matrix with (<r, 7r)- 
th element <2 a7r , where a,7T G {1, • • • , c}. Furthermore, 
define := A(x) —A(x), where A is the diagonal matrix 
whose (p, p)-th element is equal to the sum of the ele- 
ments of the p-th column of A. Hereafter we refer to L(x) 
as the weighted Laplacian of the complex graph associated 
with the given network with species concentration vector 
x. It can be verified that the sum of the elements of each 
column of L is equal to zero and that the vector Bv(x) is 
equal to —L(x)Exp(Z T Ln(x)). From equation (1), it fol- 
lows that the dynamics of enzymatic chemical reaction 
networks can be compactly written as 



x = —ZL(x)Exp (z T Ln(x)^ + Zv b {x) 



(8) 



A similar expression of the dynamics corresponding to 
mass action kinetics, in less explicit form, was obtained 
in [31]. Note that although in equation (8), Exp(Z r Ln(#)) 
is not defined if some of the components of x are equal 
to zeros, the limit at such an x exists and is used in 
equation (8) whenever some of the components of x are 
equal to zero. 

Model reduction 

For many purposes one may wish to reduce the number 
of dynamical equations of a chemical reaction network in 
such a way that the behaviour of a number of key metabo- 
lites is approximated in a satisfactory way. We propose a 
novel method for model reduction of chemical reaction 
networks governed by enzyme kinetics. Our method is 
inspired by the Kron reduction method for model reduc- 
tion of resistive electrical networks described in [25]; see 
also [32]. 

Description of the method 

Consider again a reaction network with boundary fluxes 
and with dynamics modelled by equation (8). Our model 
reduction method is based on reduction of the complex 
graph associated with the network. Let V denote the set 



of vertices of the complex graph. Reduction of the model 
is carried out by deleting certain complexes in the complex 
graph, resulting in a reduced complex graph. Deletion of a 
complex is equivalent to imposing the complex balancing 
condition on it, i.e., the condition that the net inflow into 
the complex is equal to the net outflow from it. Consider 
a subset V 0 C V of dimension c — c that we wish to delete 
in order to reduce the model. Consider a partition of L(x) 
given by 



L(x) = 



In (*) Ln(x) 
L 2 i(x) L 22 (x) 



where L u (x) G R £x£ , L 12 (x) G W x ^ c \ L 21 (x) G 
R (c-c)xc j L22{pc) e R (c-c)x(c-c) and Vq corresponds to the 

last part of the indices (denoted by 2). Consider corre- 
sponding partitions of Z and v b given by 



Z=[Zi Z 2 ]; v b (x) 
in order to rewrite (8) as 



[v bl (x) v b2 (x)f 



v bl (x) 



-Z 



L n (x) L 12 (x) 
L 2 i(x) L 22 (x) 



Exp(Zfhn(x)) 
Exp(ZjLn(#)) 



Define P :=[I^ — ^12^22 1 an( ^ ^(*) denote the Schur 
complement of L(x) with respect to the indices corre- 
sponding to V 0 . Consider now the auxiliary dynamical 
system 







v bl (x) 













L u (x) L u (x) 
L 21 (x) L 22 (x) 



W\ 

w 2 



Note that complex balancing condition on the complexes 
in V 0 can be imposed by setting the constraint y 2 = 0. This 
results in the equation 

w 2 = -L 22 (x)~ l (v b2 (x) -L 2 i(x)wi), 

leading to the reduced auxiliary dynamics 

y x — Pv b (x) — L(x)W\. 

Substituting w\ = Exp(Zf Ln(#)) in the above equation 
and making use of x = Z\ y\ -\-Z 2 y 2 = Ziyi, we then obtain 
the reduced model given by 



x = Z\ (Pv b (x) — L(#)Exp(zf Ln(^))^ . 



(9) 



Note that the reduced model is independent of the order 
of deletion of complexes. From the following Proposition, 
it follows that L(x) satisfies all the properties of a weighted 
Laplacian matrix of a reaction network corresponding to 
a complex graph with vertex set V — V 0 . 

Proposition L Consider a chemical reaction network 
with weighted Laplacian matrix L(x) ef xc correspond- 
ing to the concentration vector x. Let V denote the set of 
vertices of the complex graph associated with the network. 
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Then for any subset of vertices V r C V, the Schur comple- 
ment L(x) ofL(x) with respect to the indices corresponding 
to V r satisfies the following properties: 

1. All diagonal elements ofL(x) are positive and 
off-diagonal elements are nonnegative forx e R+. 

2. lTL(x) = 0, where c := c - dim(V r ). 

Proof. (1.) Follows from the proof of ([33], Theorem 
3.11); see also [32] for the case of symmetric L. (2.) With- 
out loss of generality, assume that the last c — c rows and 
columns of L(x) correspond to the vertex set V r . Consider 
a partition of L(x) given by 

£(*) = [w*w 12 !11 

[_L 2 i(x) L 22 (x) _ 

where L n (x) e Wf c ,L 12 (x) e W x(c ~ c \ L 21 (x) e R( c ~ c ^ c 
and L 22 (x) e r(c-c)x(c-c) # It is easy to see t j iat 

L(x) = L n (x) -Ii2W^22W _1 ^2i(^) 
Since 1 T c L{x) = 0, we obtain 

lTL n (x) + l c r _ a I 2 i (x) = t?L 12 (x) + \ T c _ l L 22 (x) = 0 
Using the above equations, we get 
tTL(x) = tT(L 11 (x)-L 12 (x)L 22 (x)- 1 L2i(x)) 

= ~^c-M (*) + l c r -a^22(^22W"^21 (X) = 0 

□ 

This proves that (9) describes the dynamics of a chemi- 
cal reaction network governed by enzyme kinetics, with a 
reduced number of complexes and with, in general, a dif- 
ferent set of boundary fluxes and reactions (edges of the 
complex graph). An appropriate choice of V 0 will ensure 
that some of the elements of x have derivative zero in (9) 
leading to a reduced number of state variables. 

The principle behind our model reduction method is to 
couple the dynamics of certain complexes to the dynam- 
ics of the neighbouring complexes or the complexes with 
which they are connected by reactions. This is done by 
complex balancing or by equating the net rate of inflow 
into the complex to the net rate of outflow from it. Note 
that for a reasonably good model reduction, it is important 
to make the right choice of complexes to be deleted. In 
the next section, we describe a procedure for making the 
choice of complexes to be deleted. Below, we illustrate our 
complex balancing procedure for model reduction with 
the help of an example. 

Example 1. We consider an example of a simple reversible 
reaction network and apply the model reduction proce- 
dure described in the paper This reaction network is shown 
below: 



For i = 1, . . . , 6, let X[ denote the concentration ofX[. For 
i = 1, ... ,4, let K]£ denote the Michaelis constant of Xi 
for the first reversible reaction. For i = 3, ... ,6, let K^ 1 
denote the Michaelis constant ofXifor the second reversible 
reaction. Let kj and k l r denote the forward and reverse rate 

constants of the first reaction and let I<j and k 1 / denote 
the respective rate constants of the second reaction. Note 

T V f T 

that k) := u 7 I>2 where Vi denotes the maximum rate 

J Km K m J 

of the first reversible reaction in the forward direction and 
the other rate constants are similarly defined. Define x := 
[x\ x 2 ... x 6 ] T , 



pi(x) 
Pi(x) 



4- 



X\ 

1,1 



X3 



+ 



X3 
1,3 



11,3 



+ 



Ki 

X5 
■11,5 



1 + 



X 2 
1,2 



#4 



+ 



#4 
1,4 



KL 

X6 



Ki 



-11,4 



Kr 



11,6 



As described for the example in the section "General 
kinetics", possible rate equations for the network are given 



by vif{x) = — 

k 11 : 



KXlX \ Vlr (x) = 



pi(x y rir\"J - pi{x y V 2 f{X) - p2(x) 

v 2r {x) = k p*(*) ' wnere v \f> v ir denote the reaction rates in 
the forward and the reverse directions respectively of the 
first reversible reaction and v 2 f, v 2r similarly denote those 
of the second reversible reaction. IfK l eq and K^ q denote the 
equilibrium constants of the first and the second reversible 
reactions respectively, then 



and 



K* =4 

eq kl 



K; 



eq 



k 11 
X 
k'J 



Now consider the reduced network 

Xi + X 2 ^ Xs + Xs 

that is obtained by deleting the complex X3 + X4 from 
the network (11). Applying the procedure described in this 
section, we first assume that the rate of inflow to the com- 
plex X3 + X4 is equal to the rate of outflow from it, i.e. 
v if — v ir = V2f — v 2r . This yields 



kjp 2 (x)xix 2 + k I r I pi(x)x 5 x 6 



X3X4, = 



kjp\(x) + k l r p 2 (x) 



(12) 



By substitution, we obtain the following expression for 
the overall rate v in the forward direction of the reduced 
network: 



X\ + X 2 ^ X3 + X4 ^ Xs + X^ 



(id 



v( x)=n f (x)-v lr ( X ) = v 2f ( X )- V2r ( X ) = q pi(x) + klpi(x) 

In the expression for p\ and p 2 in the right hand side of the 
above equation, we fix X3 and #4 at their initial values to 
obtain the following expression for v. 
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kj d %\%i — k r r ed x$X6 

^re^l ~r ~r ' T ^red,6 ' T ^red,l2 ~r T ^red,56 

Am A m A m A m A m A m 

(13) 

w/zere A™*, A?* K™ d >\K™ d > 2 , K™ d > 5 , K™ d > 6 , K™ d > 12 , K™ d > 56 
are positive constants that depend on the rate constants 
and Michaelis constants of the original network and the 
initial values ofxz and x^. We remark here that the result- 
ing equation (13) has again the form of an enzyme kinetic 
rate equation and the rates in the forward and reverse 
directions can be written in the form of equation (2) with a 
d satisfying d : -> R+. IfK r e ^ denotes the equilibrium 
constant for the reduced network, then observe that 

k l k n 

jy-red _ f f _ jA jAI 

l\y,l\y. 

For this example, our procedure yields a reduced model 
with 8 parameters while the original model has 12 parame- 
ters. Moreover while the original model has 2 reactions and 
6 state variables, the reduced model has 1 reaction and 4 
state variables. 

Error integral 

We now describe an automated procedure for making 
the choice of complexes to be deleted in such a way that 
the dynamic behaviour of the reduced model is close to 
that of the original model We quantify the difference 
between the dynamical behaviour of the original and the 
corresponding reduced model under the conditions of a 
specific dynamic event by an error integral, as defined in 
this section. 

Assume that the given biochemical network is asymp- 
totically stable around a steady state. Let M.\ denote the 
set of species which we consider to be significant from the 
point of view of experimental design. This set is a sub- 
jective choice of the scientist and contains for instance 
species whose concentrations can be experimentally mea- 
sured. Moreover complexes made of species in M.\ are 
not considered for deletion in order to reduce the model. 
Let n(M\) denote the number of elements in the set 
M\. Assume that [0, T] is the time interval over which 
we are interested to observe the difference between the 
behaviours of the full and the reduced model. We define 
the error integral (/) as 




where Xi r (t) and X[f(t) denote the concentrations of the i th 
metabolite of the reduced and the full model respectively 
at time t. Observe that / is dimensionless. The value of / 
indicates the deviation of the reduced model from the full 
model averaged over time and species. 



In order to reduce a given model of a biochemical 
reaction network, we first determine the steady state con- 
centration of each of the metabolite of the network. Then 
we rank the complexes to be deleted according to the 
error integrals corresponding to the reduced model with 
the respective complex deleted and with the initial values 
of concentrations of the species of the deleted complexes 
set at their steady state values. We then make the dele- 
tion that leads to the smallest value of /. With the reduced 
model, we then repeat this procedure. Thus we follow an 
iterative procedure with each iteration consisting of two 
steps: ranking of complexes to be deleted and deletion of 
the complex which leads to the smallest error integral. We 
stop the iteration when the error integral of the reduced 
model obtained at the end of an iteration exceeds a cer- 
tain cut-off value. The reduced model obtained at the end 
of the previous iteration is then considered as the final 
reduced model. 

For the two biochemical models considered in this 
paper, namely the yeast glycolysis model and the rat-liver 
fatty acid beta-oxidation model, the cut-off value of the 
error integral for stopping the iterative model reduction 
procedure has been set at 0.1. In general, this cut-off 
value can be set according to the desired closeness of the 
reduced model to the original. 

Difference with QSSA 

The basic premise of our model reduction approach is 
similar to the one of quasi steady state approximation 
(QSSA). In our method, we assume that some of the com- 
plexes attain rapid steady states and hence we equate the 
rates of inflows to and outflows from such complexes. On 
the other hand, in the case of QSSA, it is rather some 
of the individual metabolites that are assumed to attain 
rapid steady states and hence the rates of inflows to and 
outflows from such metabolites are equated. There are 
some other subtle differences between our approach and 
QSSA even for the case when the complexes are given by 
the species, i.e., Z is an identity matrix, as shown in the 
following example. Consider the reaction network 

X\ ^ X2 ^ X3 

with reactions governed by reversible Michaelis -Menten 
kinetics. As earlier, for i = 1, 2, 3, let X{ denote the con- 
centration of X{. For i = 1, 2, let denote the Michaelis 
constant of X{ for the first reversible reaction. For i = 2, 3, 
let K^ 1 denote the Michaelis constant of X/ for the second 
reversible reaction. Let kj and k\ denote the forward and 

reverse rate constants of the first reaction and let kj and 

I<l l denote the respective rate constants of the second reac- 

t v} T 
tion. Note that k f := -^fr where VI denotes the maximum 

J Km J 
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rate of the first reversible reaction in the forward direction 
and the other rate constants are similarly defined. Define 



Pl(xi,x 2 ) := (l + ^j + ^j 



X 2 X 3 

P2(x 2 ,x 3 ) := 1 + — ^ + 



K n,2 1 K u,3 / • 

Now assume that X 2 attains rapid steady state which, as in 
the previous example, implies that 



X 2 



kjp 2 (x 2 ,x 3 )xi +kl l pi (xi,x 2 )x 3 
kjpi(xi,x 2 ) + k l r p 2 (x 2 ,x 3 ) 



If we use QSSA method, we first need to solve for x 2 from 
the resulting quadratic equation given by 



4 [A • e>) ••»••(*" ' 



-"> 1 ( 1+ ^)-^ 3 ( 1+ 3 r ) =0 ' 



Let / be a function of two variables such that x 2 = f(xi f x 3 ) 
defines an admissible (real and positive) solution of the 
above equation. Then the resulting reaction rate in the 
forward direction (v) of the reduced network 



X\ ^ X 3 
is given by 

v(x) = 



(15) 



kjkfxi - k l r kj l x 3 



kjpi(xi,f(xi,x 3 )) + k l r p 2 (f(xi,x 3 ),x 3 ) 

The kinetics of the reduced network obtained by applica- 
tion of QSSA in this case is no longer enzyme kinetics. In 
contrast, our model reduction method yields the simple 
expression 



v(x) = 



kf d xi - k r r ed x 3 

1 _|_ *1 _|_ x 3 

L ~T ^red,l ~r ^red,3 
A m Am 



as the overall reaction rate in the forward direction of the 
reduced network (15). 
If we now consider the following reaction network 

Xi + X 2 ^ 2X3 ^ X4 + X5 

governed by the same kind of enzyme-kinetic rate laws 
as in Example 1 and assume that X 3 attains rapid equilib- 
rium, application of QSSA is more complicated as com- 
pared to the previous case as it involves solution of a 
quartic equation. On the other hand, our method leads 
to a simple expression for the overall reaction rate in the 
forward direction (v) of the reduced network 



v(x) 



kjr ed X\%2 — k r r ed X4Xs 



1 _L *1 _L x 2 1 *4 I *5 1 x l x 2 1 *4*5 ' 
± "T" r^red,l * Tr red,2 * r^red,4 T ^red.5 * ^red.12 T ^red.45 
J\m -<Mn -<Mn -<Mn 



When the complex to be deleted in a reaction net- 
work is made up of more than one species, QSSA cannot 



be applied in some cases. Example 1 is one such case. 
With reference to this example, using equation (12), one 
needs to solve for x 3 and #4 in order to apply QSSA. 
However, this is not possible as there is one equation 
and two unknowns x 3 and #4 to be solved for. Thus our 
model reduction method is more effective in dealing with 
deletion of complexes made up of more than one species. 

Effect of model reduction 

In this section, we show the graph restructuring of our 
reduced model in terms of its corresponding full model 
for three particular types of linkage classes of biochemi- 
cal reaction networks. Note that the deletion of a set of 
complexes in one linkage class does not affect the mathe- 
matical models of the reactions of the other linkage classes 
of the network. 
Type 1 linkage class: 



Full Network: 



C\ ^ C 2 ^ C 3 ^ 



(16) 



Reduced Network: C\ ^ C 3 ^ ^ C n (17) 

Consider a linkage class with reversible reactions occur- 
ing between consecutive elements of the set of distinct 
complexes {Ci, C 2 , . . . , C n } as in (16). The reduced network 
obtained by deleting the complex C 2 is given by (17), where 
the two reversible reactions, C\ ^ C 2 and C 2 ^ C 3 of the 
full network are replaced by a reversible reaction C\ ^ C 3 
in the reduced network. If the kinetics of the reactions 
C\ ^ C 2 and C 2 ^ C 3 are of the same type, the reaction 
C\ ^ C 3 of the reduced network has the same kinetics. 
This has been shown for a particular type of reactions in 
Example 1. The rate equations of all the remaining reac- 
tions of the reduced network are the same as those of the 
full network. 

A special case of linkage class (16) is the following: 



(18) 



Deletion of the complex C 2 in this case is equivalent to 
deletion of the linkage class from the network. Such a 
deletion provides a close approximation to the original 
network if the reaction (18) has very little effect on the 
dynamics of the network. 
Type 2 linkage class: 

• • • C\ ^ c 2 -> c 3 



Full Network 1: 



Reduced Network 1: 



Full Network 2: 



C\ -> c 2 ^ c 3 



Reduced Network 2: • • • C\ -> C 3 



(19) 



(20) 



(21) 



(22) 
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In this type of a linkage class, we have one complex C2 
involved in both a reversible reaction and an irreversible 
reaction as shown in (19) and (21). The reversible and the 
irreversible reactions in which C2 is involved in the full 
network are replaced with a single irreversible reaction 
in the reduced network as in (20) and (22). The deletion 
of C2 does not affect the mathematical description of the 
remaining reactions of the linkage class. As with the pre- 
vious type of linkage class, if the kinetics of the reversible 
and irreversible reactions of the full network are of the 
same type, the reaction C\ C3 of the reduced network 
has the same kinetics. 
Type 3 linkage class: 

Full Network 



Results and discussion 

Yeast Glycolysis model 

We have applied our model reduction procedure on a 
detailed kinetic model of yeast glycolysis [26] . A schematic 
representation of the model is shown in the left-hand 
panel of Figure 2. The corresponding detailed mathemat- 
ical model can be found in [26] . This network consists of 
type 1 and type 2 linkage classes. 

For the modelling, we have considered the following as 
external fluxes as indicated in Figure 2: 

1. uptake of extracellular glucose (Glco) into the cell; 

2. conversion of trehalose into intracellular glucose 
(Glci); 

3. production of trehalose from glucose 6-phosphate 
(G6P); 

4. production of glycerol from TRIO (TRIO is a pool 
summing up dihydroxyacetone phosphate and 
glyceraldehyde 3-phosphate); 

5. production of succinate from pyruvate (PYR); 

6. production of acetate from acetaldehyde (Ac Aid); 

7. production of ethanol from AcAld. 



Reduced Network 

- d 4 c 3 - 

C4 

i (24) 

In this type of linkage classes, we have a complex that is 
involved in more than two reactions as the one shown in 
(23). In this example, if we delete the complex C2, we arrive 
at the reduced network shown in (24). This deletion has no 
effect on the mathematical description of the remaining 
part of the network. Similar to the previous type of linkage 
classes, the reduced model of this type of a linkage class 
inherits the kinetics of the full model. 

Notice that though the number of complexes in the 
reduced network is less than that of the original network, 
the number of reactions is the same in both. The reaction 
constants of reaction 4 of the reduced network depends 
on the reaction constants of reactions 1 and 2 of the full 
network; those of reaction 5 depend on those of reactions 
1 and 3 and those of reaction 6 depend on those of reac- 
tions 2 and 3. Furthermore, deletion of the complex C2 in 
this case does not lead to a reduction in the number of 
parameters. 

In the next section, we show the application of our 
model reduction method on a yeast glycolysis model and 
a beta oxidation model in rat liver. In both these cases, 
we encounter linkage classes belonging to one of the three 
types mentioned above. 



The model reduction procedure is applied to a glucose 
upshift as described in [26] and the corresponding param- 
eter set was chosen [26]. Under these conditions, the 
imposed concentrations of Glco are 0.2 mM and 5 mM for 
t < 0 and t > 0 respectively, and as observed experimen- 
tally the corresponding concentrations of ATP are equal 
to 5 mM and 2.5 mM. It is assumed that the network is at 
steady state for t < 0. It is found that the model is asymp- 
totically stable. The reactions of the network are governed 
by enzyme kinetics and we can write the equations of 
the model in the same form as equation (8). The set of 
significant species (.Mi) for the model consists of Glci, 
TRIO, BPG, PYR, AcAld and NADH. Figure 3 depicts the 
minimum error integrals as a function of the number of 
deleted complexes. Since deletion of a sixth complex leads 
to the error integral exceeding 0.1 in value, we stop the 
model reduction process after deletion of five complexes. 
The order of deletion is F6P, G6P, 2-phosphoglycerate 
(P2G), 3-phosphoglycerate (P3G) and phosphoenoylpyru- 
vate (PEP). Deletion of P3G and P2G produces the same 
effect as the deletion of complex C2 from a type 1 link- 
age class as described in the section on model reduction. 
Deletion of each of the remaining complexes produces the 
same effect as the deletion of complex C2 from a type 2 
linkage class. Each deletion results in the shortening of 
the chain length of the linkage class to which the complex 
belongs (see Figure 2, right-hand panel). 

It is observed that there is a good agreement between 
the transient behaviours of most of the metabolites when 
comparing the original model to the reduced model with 
five complexes deleted. The main reason for this is that 
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Figure 2 Schematic of the original and reduced yeast-glycolysis networks. The left-hand panel is a schematic representation of the yeast 
glycolysis model used for model reduction. The full model description and an explanation of all the abbreviations is found in [26]. The right-hand 
panel represents the reduced model after deleting 5 complexes (F6P, P2G, P3G, G6P and PEP). The blue arrows in the two panels indicate external 
fluxes. 
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Figure 3 Reduction of the yeast glycolysis model. Left-hand panel: minimum error integral versus number of deleted complexes (we have taken 
T = 1 .5 min for the computation of the error integral). Right-hand panel: identity of the deleted complexes and their convergence times. 
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the reactions of the original model that are missing in the 
reduced one (see Figure 2) are the close-to-equilibrium 
reactions of the network. The concentration of F16BP has 
a strong effect on enzyme rates, both locally on the reac- 
tions in which it is involved, and distantly on the rate 
of pyruvate kinase (PYK). Together, this provides a clear 
explanation why F16BP ends up last in the order of com- 
plexes to be deleted (see Figure 3). The table of Figure 3 
gives the convergence times of the six metabolites (com- 
plexes) that were considered for deletion. These are the 
times that it takes for the metabolite concentrations to 
achieve 95% of their concentration change going from one 
steady state to the other. Observe that the order of deletion 
of the complexes is not the same as the increasing order of 
their convergence times. 

While the original model has 12 variables, 88 parame- 
ters and 12 reactions, the reduced model has 7 variables, 
50 parameters and 7 reactions. We have provided the Mat- 
lab files corresponding to the original model, the reduced 
model and automation of the model reduction proce- 
dure for the yeast glycolysis model, as additional files 
(see Additional files 2, 3 and 4 respectively). Additional 
file 4 gives as output the order of deletion of com- 
plexes and the error integral at each step of deletion. 
Additional files 5, 6, 7 and 8 are Matlab files that are 
required to run Additional file 4. We have submitted to 



Biomodels two SBML files corresponding to the origi- 
nal and the reduced yeast glycolysis models (submission 
identifiers MODEL1403250001 and MODEL1403250002 
respectively). The evolution of the concentrations of 
Glci, PYR, TRIO, ACALD, BPG and NADH is com- 
pared between the different stages of model reduction in 
Figure 4. 

Fatty acid beta-oxidation model 

Subsequently, we have applied our model reduction 
method to a very different type of biochemical network, 
i.e., a model of fatty-acid beta-oxidation in rat liver [27]. 
This network mostly consists of type 1 and type 3 link- 
age classes. The model is shown schematically in Figure 5. 
Acyl carnitines of even chain lengths are broken down in 
the mitochondria to produce acetyl CoA. The acyl car- 
nitines are first converted to acyl Co As of the same chain 
length within the mitochondria by the enzyme carnitine 
palmitoyl transferase II. Through a series of enzymatic 
reactions involving the enzymes acyl CoA dehydrogenase, 
crotonase (CROT), hydroxyacyl CoA dehydrogenase 
(M/SCHAD), ^-ketothiolase (MCKAT) and mitochon- 
drial trifunctional protein (MTP), the carbon chain 
length of each molecule of acyl CoA is shortened by 
2 at the same time producing one molecule of acetyl 
CoA. 
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Carnitine CoASH 





Figure 5 Schematic of a model of beta-oxidation in rat liver. The ellipses within, on and outside the boundary of the mitochondrion represent 
the enzymes responsible for the beta-oxidation. The full model description and explanation of all the abbreviation can be found in [27]. 



The carbon chain-lengths of the compounds processed 
by each of these enzymes is indicated within the ellipses. 
The full model has 42 state variables, 160 parameters and 
56 reactions. The mathematical model is first written in 
the form (8). It is found that the model is asymptoti- 
cally stable. At steady state, each enzyme has a constant 
reaction flux (rate). 

By making use of the biochemical knowledge of the beta 
oxidation model, an initial reduction of the beta oxida- 
tion model has been performed as described below. This 
procedure holds only for this particular model and is not 
obtained by an automated deletion of complexes. It can be 
seen from Figure 5 that certain acyl Co As having the same 
carbon chain lengths can be dehydrogenated by multiple 
enzymes, for example C8 acyl CoA can be dehydrogenated 
by three enzymes MCAD, LCAD or VLCAD. By checking 
the steady state fluxes of the same reaction through these 
different enzymes, we can reduce the model. For exam- 
ple, we found that the flux of dehydrogenation of C6 acyl 
CoA by MCAD is 99 times that by SCAD. Thus the model 
can be reduced by assuming that C6 acyl CoA is dehydro- 
genated by MCAD alone and not by SCAD. Similarly, we 
found that the model can be further reduced by assum- 
ing that C4 acyl CoA is dehydrogenated by SCAD alone 
and C8 acyl CoA is dehydrogenated by MCAD alone as 
the dehydrogenation fluxes of these compounds via other 
enzymes is comparatively very low. 



Subsequently, we performed a further reduction of the 
beta oxidation model by deletion of certain complexes 
according to the procedure described earlier. It can be 
seen from Figure 5 that the shortening of chain length of 
acyl CoAs of chain lengths 16, 14, 12, 10 and 8 can happen 
via two routes, i.e. either via the enzyme MTP or via the 
sequence of enzymes crotonase, M/SCHAD and MCKAT. 
Deletion of the complexes hydroxyacyl and ketoacyl Co As 
of even chain lengths between 8 and 16 is equivalent to 
having all acyl Co As of chain lengths 8 till 16 reduced 
by only the first route. It is found that such a deletion 
produces a reduced model which has a very similar tran- 
sient behaviour as the full model. The observation that the 
steady state fluxes through the first route is about a hun- 
dred times larger than the steady state fluxes through the 
second route explains the fact that removal of the second 
route does not have much effect on the dynamics of the 
rest of the system. 

Note that the shortening of chain length of acyl CoAs 
of chain lengths 6 and 4 happens only via the sequence 
of enzymes crotonase, M/SCHAD and MCKAT. It is 
found that the combined effect of enzymes crotonase and 
M/SCHAD can be produced by a single fictitious enzyme 
which we have termed as CRMS. This is done by deleting 
the complexes C4 hydroxyacyl CoA and C6 hydroxyacyl 
CoA. Again such a deletion has a negligible effect on the 
dynamics of the rest of the system. 
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Each of the 12 deletions of complexes mentioned above 
produces the same effect as the deletion of the complex 
C2 from a type 1 linkage class as described in the section 
on effect of model reduction. Deletion of hydroxyacyl and 
ketoacyl CoAs of even chain lengths between 16 and 8 is 
equivalent to removal of the linkage classes to which they 
belong. Deletion of C4 and C6 hydroxyacyl CoA is equiv- 
alent to reducing the number of reactions in the linkage 
class to which they belong, which in turn is equivalent to 
replacing the enzymes crotonase and M/SCHAD with a 
single enzyme CRMS. 

The set of significant species (.Mi) for the model con- 
sists of all the acyl carnitines in the cytoplasm and the 
mitochondria, since these have been measured for the 
original model validation [27] . Figure 6 depicts the values 
of minimum error integrals vs the number of complexes 
deleted in order to obtain a reduced model of beta oxi- 
dation. It can be seen that when the fourteenth complex 
is deleted from the model, the minimum error integral 
is more than 0.1. Therefore as discussed earlier, we stop 
our model reduction procedure after deletion of 13 com- 
plexes. The first 13 complexes that are deleted are the ones 
that were mentioned in the previous paragraphs in addi- 
tion to C8 acyl CoA. Deletion of C8 acyl CoA produces the 
same effect as the deletion of the complex C2 from a type 
3 linkage class as described earlier. This deletion does not 
lead to a reduction in the number of reactions or param- 
eters; only the number of complexes is reduced. In the 
following we show the local effect of deletion of C8 acyl 
CoA schematically. 
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The reduced model obtained by incorporating all the 
deletions mentioned above except the deletion of the com- 
plex C8 acyl CoA is depicted in Figure 7. It has 31 state 
variables, 118 parameters and 37 reactions while the full 
model has 43 state variables, 160 parameters and 56 reac- 
tions. It is found that the transient behaviour of all the 
state variables of the reduced model with 12 complexes 
deleted are in good agreement with those obtained using 
the full model. Both the original and the reduced model 
have been provided as additional files (see Additional 
files 9 and 10 respectively). An SBML file corresponding 
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Figure 6 Minimum error integral vs number of complexes deleted from the beta-oxidation model. We have taken T = 25 min for the 

computation of the error integral. 
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to the reduced model has been submitted to Biomod- 
els (submission identifier MODEL1403250000). Figure 8 
depicts the comparison of the transient behaviours of the 
concentrations of all the acyl carnitines in the cytoplasm. 

Discussion 

Our model reduction method involves rewriting of the 
complex graph corresponding to a given network. In the 
literature, there are other model reduction methods that 
involve graph rewriting like QSSA and quasi-equilibrium 



procedures, the method of [21] and the graph rewriting of 
monomolecular reaction networks described in [3]. The 
difference between these methods and ours is that these 
involve rewriting of the species-reactions graph unlike 
ours where the complex graph is rewritten. 

In order to come up with a reasonably good reduced 
model, we follow an iterative procedure as described 
in the subsection titled "Error integral". This iterative 
procedure involves computation of error integrals of a 
number of reduced models which are obtained by deletion 
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Figure 8 Comparison of concentration profiles of acyl carnitines between the full and reduced beta-oxidation models. 
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of different combinations of complexes. In the beta- 
oxidation example, where the original model consists of 43 
species and the reduced one has 30 species, the required 
number of simulations is (43 —M) + (42 - M) + (41 —M) + 

h (31 —M) where M is the number of species that must 

still be present in the reduced model (and hence, are not 
considered in the subset of complexes for deletion). In our 
example, M = 14 (corresponding to all Acyl-carnitines in 
the cytosol and mitochondria) and the total simulations 
using our procedure are 299. More precisely, if there are 
N number of species, M number of important species and 
if L is the dimension of the reduced species, then the total 
number of simulations is ^ L (/-M) = -M)(N- 
L + 1). Although our procedure is time consuming, it is 
simple and systematic, ensures that the reduced model 
closely mimics the transients of the original model and is 
scalable to a large model. We are currently investigating 
the computation of error bounds based on the algebraic 
property of the Laplacian matrices and we foresee that this 
knowledge might allow us to directly obtain the optimal 
set of complexes to be deleted. 

The cut-off value of the error integral at which we 
stop our iterative model reduction procedure can be set 
according to the desired closeness to the original model. 
For the two models discussed in this paper, this value has 
been set at 0.1 and it has been found that if this value is set 
at 0.2, then the resulting comparison plots of the transient 
behaviours of the original vs the reduced models reveal 
significant visual differences. Another method of coming 
up with a cut-off value for the error integral is to look for 
a sudden, steep increase in the error integral histogram. 
However in the case where the error-integral histogram 
shows uniform increase with an increase in the number of 
complexes deleted, this method cannot be used. 

Note that in principle, complexes can be deleted to 
reduce any biochemical model, whether asymptotically 
stable or not. However, for the computation of error inte- 
gral as described in the paper, it is necessary that the 
original network is asymptotically stable around a steady 
state, since the computation makes use of the steady state 
concentrations of some species of the network. 

In some cases, deletion of certain complexes can drasti- 
cally change the dynamics of the system. An example is the 
following. Assume thatXi,X2, . . . ,Xio are distinct species 
of the network. Consider the following reaction network 
consisting of three reversible reactions: 

X\ + X2 ^ X3 + X4, 
X4 + Xs ^ X$ + X] 
X7 + X$ ^ X9 + X\o 

Now consider deletion of the complex X4 + X5 from 
the above reaction network. This deletion leads to the 



deletion of the second reaction of the network. Since the 
first and the third reaction of the network do not have 
common species, deletion of the second reaction leads to 
the first and the third reactions occuring independently of 
each other. This will lead to the reduced model exhibit- 
ing a behaviour which is not close to the behaviour of the 
original model. Although one can quantify the difference 
between the behaviours using error integrals, avoiding 
such deletions will reduce the computational effort in 
making the choice of complexes to be deleted. Examples 
of such deletions are deletion of complexes TRIO or BPG 
in the yeast glycolysis model and deletion of complex C16 
Acyl Carnitine in the beta oxidation model. 

We remark here that one should update a reduced 
model when major changes are made to the original 
model. For example, if we consider the rat liver beta- 
oxidation model of a MCAD-deficient animal, this model 
will not have the enzyme MCAD. In this case, the dynamic 
behaviour of the reduced model of Figure 7 may be far 
from that of the original model of Figure 5 without the 
enzyme MCAD. 

Conclusions 

In this paper, we have outlined a method for model reduc- 
tion of biochemical reaction networks that are asymptot- 
ically stable around a steady state. The principle behind 
our model reduction method is to couple the dynamics 
of certain complexes to the dynamics of the neighbouring 
complexes in such a way that the net rate of inflow into 
the complex is equal to the net rate of outflow from it. 
We provide an algorithm for choosing the complexes to be 
deleted in such a way that the transient behaviour of the 
significant species of the reduced model is close to that of 
the original model. Apart from a reduction in the number 
of state variables, our model reduction method also leads 
to a reduction in the number of reactions and parameters 
of the model. This in turn leads to an improved com- 
putational effort required in order to analyze the model. 
Our model reduction procedure ensures that the reduced 
model mimics the full model well under the conditions of 
a specific dynamic event. A different reduced model with 
a different set of deleted complexes may be produced by 
our procedure under the conditions of a different dynamic 
event. 

In Additional file 1, we give a list of some well-known 
enzyme kinetic rate laws for which our model reduction 
method is applicable. Some of the rate laws provided in 
this list are rate laws for irreversible reactions. However 
our method is also applicable for the reversible version 
of these rate laws. This is because a reversible reac- 
tion can be split into its forward and reverse reaction 
components and if our method is applicable to each 
of these components then it is applicable to the given 
reversible reaction. For example, since our method is 
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applicable to reactions governed by irreversible mass 
action kinetics, it is also applicable for reactions gov- 
erned by reversible mass action kinetics. The list pro- 
vided in Additional file 1 is not an exhaustive list of all 
enzyme kinetic rate laws for which our method is appli- 
cable. We have checked that our method is applicable 
to all the rate laws that are included in the COPASI 
software [34]. 

We have applied our method on a yeast glycolysis model 
and a beta-oxidation model in rat liver and have observed 
a good agreement between the transient behaviours of 
the original and the reduced models. In these cases, our 
model reduction method has also helped us significantly 
in improving our understanding of the dynamics of the 
networks, for example by identification of the close-to- 
equilibrium reactions of the yeast glycolysis model and 
of the reaction pathways with comparatively low fluxes 
in the rat liver beta-oxidation model. In the case of the 
yeast glycolysis model, we found that the order of deletion 
of complexes (metabolites) is not the same as the order 
of their convergence times. This shows that in this case, 
our model reduction approach will perform differently 
from a time-scale separation technique where metabolites 
with faster convergence times are assumed to attain rapid 
steady state. As observed in the two examples of metabolic 
reaction networks, our model reduction method helps in 
removing redundancies in the given network (in the yeast 
glycolysis network, redundant reactions were coupled and 
in the rat liver beta-oxidation model, redundant pathways 
were removed). 

Recently, there has been an interest in kinetic mod- 
elling of large genome-scale networks. For example, [35] 
gives a method for building a parameterized genome-scale 
kinetic model of a network consisting of 956 reactions and 
820 metabolites. We envisage that our model reduction 
method will be particularly useful to reduce the com- 
plexity of such large genome-scale kinetic models in the 
future. Likewise, it can also be useful in reducing mul- 
tiscale models of biochemical reaction networks which 
are computationally burdensome to analyze. In a large- 
scale kinetic model of a biochemical reaction network, our 
model reduction method can be used to reduce all path- 
ways except the ones of interest into which we would like 
to zoom-in and analyze. 

Currently, we are using our reduced model of the rat 
liver beta-oxidation model that we explained in an earlier 
section in order to check for possible bifurcations. In this 
respect, the reduced model provides an advantage over 
the full model as it requires lesser computational effort 
while checking for bifurcation. 

Ethics 

The authors declare that no experiments have been per- 
formed as part of the research for this manuscript. 



Additional files 



Additional file 1 : Enzyme kinetic rate laws for which our model 
reduction method is applicable. 

Additional file 2: Original yeast glycolysis model. 

Additional file 3: Reduced yeast glycolysis model. 

Additional file 4: Automation of our model reduction procedure for 
the yeast glycolysis model. 

Additional file 5: Reduced yeast glycolysis model with a prespecified 
number of deleted complexes. 

Additional file 6: A Matlab file required to run Additional file 4. 

Additional file 7: A Matlab function for computing the final steady 
state with given initial conditions. 

Additional file 8: A program to compute the error integral between 
the original and a given reduced model of the yeast glycolysis. 

Additional file 9: Original model of rat liver beta oxidation. 

Additional file 1 0: Reduced model of rat liver beta oxidation. 



Competing interests 

The authors declare that they have no competing interests. 
Authors' contributions 

AJvdS, SR and BJ conceived the idea of model reduction. AJvdS prepared a 
first draft of some of the subsections under the "Methods" section. BMB 
conceived the idea of stepwise automated reduction procedure. KvE provided 
the original models of yeast glycolysis and rat liver beta oxidation. SR and BJ 
carried out reduction of these models. SR prepared the manuscript with the 
help of BMB and BJ. All authors read and approved the final manuscript. 

Acknowledgements 

We thank Dr. Guy-Bart Stan and Juan Kuntz for their critical comments on an 
earlier draft of the manuscript. We also thank the reviewers of this manuscript 
for patiently going through it and providing suggestions to improve its 
readability. Their comments have helped us in improving the quality of the 
manuscript. This work was funded by the Netherlands Organization for 
Scientific Research (NWO) to the Systems Biology Centre for Energy 
Metabolism and Ageing. AJvdS received funding from the European Union 
Seventh Framework Programme [FP7/2007-2013] under grant agreement 
No.257462 HYCON2 Network of Excellence. KvE has been funded by the 
Netherlands Genomics Initiative via the Netherlands Centre for Systems 
Biology and by the Top Institute Food and Nutrition. BMB received a Rosalind 
Franklin Fellowship from the University of Groningen. 

Author details 

1 Systems Biology Center for Energy Metabolism and Ageing, University of 
Groningen, ERIBA, Antonius Deusinglaan 1 , 971 3 AV Groningen, Netherlands. 
2 Johann Bernoulli Institute for Mathematics and Computer Science, University 
of Groningen, P.O. Box 407, 9700 AK Groningen, Netherlands, department of 
Pediatrics, Center for Liver, Digestive and Metabolic Diseases, University 
Medical Center Groningen, University of Groningen, Hanzeplein 1, 9713 GZ, 
Groningen, Netherlands, institute of Technology and Management, 
Nijenborgh 4, University of Groningen, 9747 AG Groningen, Netherlands. 

Received: 29 October 201 3 Accepted: 23 April 201 4 
Published: 3 May 2014 

References 

1 . Chickarmane V, Paladugu SR, Bergmann F, Sauro HM: Bifurcation 
discoverly tools. Bioinformatics 2005, 21 :3688-3690. 

2. Levering J, Kummer U, Becker K, Sahle S: Glycolytic oscillations in a 
model of a lactic acid bacterium metabolism. Biophys Chem 201 3, 
172:53-60. 

3. Radulescu O, Gorban AN, Zinovyev A, Noel V: Reduction of dynamical 
biochemical reaction networks in computational biology. Front 
Genet 2012,3:00131. 



Rao et al. BMC Systems Biology 201 4, 8:52 Page 1 7 of 1 7 

http://www.biomedcentral.eom/1752-0509/8/52 



10. 



12. 



13. 



14. 



15. 



16. 



17. 



19. 



20. 



21. 



22. 



23. 



24. 



25. 
26. 



27. 



28. 



29. 



Radulescu 0, Gorban AN, Zinovyev A, Lilienbaum A: Robust 
simplifications of multiscale biochemical networks. BMC Syst Biol 
2008, 2:86. 

Roussel MR, Fraser SJ: Invariant manifold methods for metabolic 

model reduction. Chaos 2001, 1 1:196-206. 

Gorban AN, Karlin IV: Method of invariant manifold for chemical 

kinetics. Chem Eng Sci 2003, 58:4751-4768. 

Sunnaker M, Cedersund G, Jirstrand M: A method for zooming of 

nonlinear models of biochemical systems. BMC Syst Biol 201 1 , 5:1 40. 

Nikerel IE, van Winden WA, Verheijen PJT, Heijnen JJ: Model reduction 

and a priori, kinetic parameter identifiability analysis using 

metabolome time series for metabolic reaction networks with linlog 

kinetics. Metab Eng 2009, 1 1 :20-30. 

Segel IH: Enzymes. In Biochemical Calculations: How to Solve 

Mathematical Problems in General Biochemistry, 2nd edition. New York: 

Wiley; 1968:208-323. 

Segel LA, Slemrod M: The quasi-steady-state assumption: a case 

study in perturbation. SIAM Rev SoclndAppI Math 1989, 31:446-477. 
Hardin HM: Handling biological complexity: as simple as possible but 
not simpler. Ph.D. Thesis. Vrije Universiteit, Amsterdam, 2010. 
Surovtsova I, Sahle S, Pahle J, Kummer U: Approaches to complexity 
reduction in a systems biology research environment (SYCAMORE). 
In Proceedings of the 2006 Winter Simulation Conference, 3-6 December 
2006. Edited by Perrone LF, Wieland FP, Liu J, Lawson BG, Nicol DM, 
Fujimoto RM. Monterey, California, USA: IEEE; 2006:1683-1689. 
Kourdis PD, Palasantza AG, Goussis DA: Algorithmic asymptotic analysis 
of the NF-kB signaling system. ComputMath Appl 201 3, 65:1 5 1 6-1 534. 
Noel V, Grigoreiv D, Vakulenko S, Radulescu O: Tropical geometries and 
dynamics of biochemical networks. Application to hybrid cell cycle 
models. Electron Notes TheorComputSci 201 2, 284:75-91 . 
Liu G, Swihart MT, Neelamegham S: Sensitivity, principle component 
and flux analysis applied to signal transduction: the case of 
epidermal growth factor mediated signaling. Bioinformatics 2005, 
21:1194-1202. 

Maurya MR, Bornheimer SJ, Venkatasubramanian V, Subramaniam S: 
Reduced-order modelling of biochemical networks: application to 
the GTPase-cycle signalling module. IET Syst Biol 2005, 152:229-242. 
Apri M, de Gee M, Molenaar J: Complexity reduction preserving 
dynamical behaviour of biochemical networks. J TheorBiol 2012, 
304:16-26. 

Androulakis IP: Kinetic mechanism reduction based on an integer 

programming approach. AlChE J 2000, 46:361 -371 . 

Bhattacharjee B, Schwer DA, Barton PI, Green WH: Optimally reduced 

kinetic models: reaction elimination in large-scale kinetic 

mechanisms. Combust Flame 2003, 1 35:1 9 1 -208. 

Petzold L, Zhu W: Model reduction for chemical kinetics: an 

optimization approach. AlChE J 1 999, 45:869-886. 

Clarke BL: General method for simplifying chemical networks while 

preserving overall stoichiometry in reduced mechanisms. J Chem 

Phys] 992, 97:4066-4071. 

Dano S, Madsen MF, Schmidt H, Sedursund G: Reduction of a 
biochemical model with preservation of its basic dynamics 
properties. FEBSJ 2006, 273:4862-4877. 

Schmidt H, Madsen MF, Dano S, Sedursund G: Complexity reduction of 
biochemical rate expressions. Bioinformatics 2008, 24:848-854. 

Bollobas B: Modern Graph, Theory. Graduate Texts in Mathematics 184. New 
York: Springer; 1998. 

Kron G: Tensor Analysis of Networks. New York: Wiley; 1 939. 
van Eunen K, Kiewiet JAL, Westerhoff HV, Bakker BM: Testing 
biochemistry revisited: how in vivo metabolism can be understood 
from in vitro enzyme kinetics. PLoS Comput Biol 2012, 8(4):e 1002483. 
van Eunen K, Simons SM, Gerding A, Bleeker A, den Besten G,TouwCM, 
Houten SM, Groen BK, Krab K, Reijngoud DJ, Bakker BM: Biochemical 
competition makes fatty-acid /? -oxidation vulnerable to substrate 
overload. PLoS Comput Biol 201 3, 9(8):e1 003 1 86. 
Horn FJM, Jackson R: General mass action kinetics. Arch Rational Mech 
Anal 1972, 47:81-116. 

Horn FJM: Necessary and sufficient conditions for complex balancing 
in chemical kinetics. Arch Rational Mech Anal 19/ 7 2, 49:172-186. 



30. Feinberg M: Chemical reaction network structure and the stability of 
complex isothermal reactors -I. The deficiency zero and deficiency 
one theorems. Chem Eng Sci 1 987, 43:2229-2268. 

3 1 . Sontag ED: Structure and stability of certain chemical networks and 
applications to the kinetic proofreading model of T-cell receptor 
signal transduction. IEEE Trans Autom Control 2001 , 46:1 028-1 047. 

32. van der Schaft, A J: Characterization and partial synthesis of the 
behavior of resistive circuits at their terminals. Syst Contr Lett 201 0, 
59:423-428. 

33. Niezink NMD: Consensus in networked multi-agent systems. Master's 
thesis in Applied Mathematics. Faculty of Mathematics and Natural 
Sciences, University of Groningen; 201 1. 

34. Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, Singhal M, Xu L, 
Mendes P, Kummer U: COPASI — a COmplex PAthway Simulator. 
Bioinformatics 2006, 22:3067-3074. 

35. Smallbone K, Simeonidis E, Swainston N, Mendes P: Towards a 
genome-scale kinetic model of cellular metabolism. BMC Syst Biol 
2010, 4:6. 



doi:1 0.1 186/1752-0509-8-52 

Cite this article as: Rao etal:. A model reduction method for biochemical 
reaction networks. BMC Systems Biology 2014 8:52. 



Submit your next manuscript to BioMed Central 
and take full advantage of: 

• Convenient online submission 

• Thorough peer review 

• No space constraints or color figure charges 

• Immediate publication on acceptance 

• Inclusion in PubMed, CAS, Scopus and Google Scholar 

• Research which is freely available for redistribution 



Submit your manuscript at 
www.biomedcentral.com/submit 



o 



BioMed Central 



